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ABSTRACT 

The traditional way of estimating the gravitational field from observed motions of 
test objects is based on the virial relation between their kinetic and potential energy. We 
find a more efficient method. It is based on the natural presumption that the objects are 
observed at a random moment of time and therefore have random orbital time phases. 
The proposed estimator, which we call "orbital roulette" , checks the randomness of the 
phases. The method has the following advantages: (1) It estimates accurately Kep- 
lerian (point-mass) potentials as well as non-Keplerian potentials where the unknown 
gravitating mass is distributed in space. (2) It is a complete statistical estimator: it 
checks a trial potential and accepts it or rules it out with a certain significance level; the 
best-fit measurement is thus supplemented with error bars at any confidence level. (3) 
It needs no d priori assumptions about the distribution of orbital parameters of the test 
bodies. We test our estimator with Monte-Carlo-generated motions and demonstrate 
its efficiency. Useful applications include the Galactic Center, dark- matter halo of the 
Galaxy, and clusters of stars or galaxies. 

Subject headings: galaxies: general — stars: stellar dynamics 



1. Introduction 

Estimation of a central point mass M from measured positions and velocities of its N satellites 
is a classic problem of astronomy which applies to various gravitating systems. More generally, the 
gravitating mass is distributed in space with density p(r), so that N test bodies move in an unknown 
potential 3>(r), and astronomers estimate <E> from the observed motions. In many cases the orbital 
periods of the test bodies are much longer than the age of modern astronomy and therefore only 
instantaneous positions and velocities are available rather than the full orbits of the bodies. This 
does not allow one to find the exact M even in a point-mass problem. Instead, statistical methods 
are used to obtain an approximate estimate. 

Customary methods of mass estimation stem from the virial relation between the mean po- 
tential energy and mean kinetic energy. It enables an estimation of <5 from the observed kinetic 
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energy of test bodies. One notes, however, that in a finite system of N test bodies, the instanta- 
neous relation between their potential and kinetic energies deviates from the virial relation. The 
variance of such deviations is unknown to the observer because it depends on the unknown orbital 
parameters of the bodies. For this reason, the observer does not possess a well defined error of the 
estimated mass. 

The virial theorem itself is a statistical statement about the potential and kinetic energies 
that is based on the presumption that the observed bodies have random orbital time phases. 1 Thus, 
the virial relation is not the original information we possess but a derivative of the random-phase 
principle. The latter simply states that the time we point our telescope to the system is random 
from the point of view of each body. One can make use of this basic principle directly, without 
invoking the virial theorem. In this paper, we follow this approach and develop a new method of 
potential estimation which we call orbital roulette (because the true orbital phase obeys the same 
statistics as a fair roulette, see § 3). We show that this method is principally better and in practice 
more efficient than the virial estimator. 

In many applications, the estimation of potential $ is further complicated by the lack of full 3D 
information on the positions rj and velocities Vj of test bodies: often only projections on the sky or 
line-of-sight components are measured. A projected version of orbital roulette will be addressed in a 
future paper. Here we develop the basic method in its full 3D version and discuss the astronomical 
data sets to which it is immediately applicable. 

The idea of our method is as follows. Consider a point-mass problem with Keplerian potential 
<3? = —GM/r and suppose one tries to infer the central mass M from observed motions of its 
satellites. For any trial mass M tria i, the observed positions r$ and velocities Vj of the satellites 
uniquely determine their orbits. Thus, for any M tria i, one can find the current orbital phases of 
the satellites. If M tr i a i is smaller than the true mass M true , the orbital phases will be found near 
the pericenters, and if M tr i a i > M true then the phases will cluster near the apocenters. Only for 
Atrial = ^true is t ne orbital roulette unbiased and the calculated phases are randomly distributed 
between the pericenter and the apocenter. One can therefore figure out Mt rue by checking the 
randomness of the inferred phases. In § 3, we put this idea on a quantitative basis and show how 
the lower and upper bounds on M are obtained at a given confidence level. 

In § 4 we extend the method to non-Keplerian potentials that have more than one unknown 
parameter and refine the roulette principle: the inferred orbital phases must be consistent with ran- 
dom numbers and uncorrelated with the inferred integrals of motion. This refinement is important 
for estimation of potentials created by distributed mass. 



lr The random- phase condition is satisfied in many astronomical systems including stellar clusters and clusters of 
galaxies (in planetary systems it may not hold if planets interact with each other and are locked in resonances). The 
virial relation is evidently not valid if the motions of the test bodies are correlated and the time of observation is 
chosen so that, for example, all bodies are near their pericenters. We do not consider the case of correlated test 
bodies here. 
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One problem in the previously proposed methods was the need for a priori assumptions about 
the test bodies' population. For instance, Little & Tremaine (1987) and Kochanek (1996) used 
Bayesian analysis to estimate the mass of our Galaxy. In order to make progress, they needed some 
assumptions about the Galactic satellites' distribution function. 2 The roulette method makes no 
assumption about the orbits of test bodies. For a trial 3>(r), it simply reconstructs the orbits and 
checks the obtained orbital phases. 5>(r) is deemed a good estimate at a confidence level C if the 
phases are consistent with random distribution at the confidence level C. This test extracts the 
maximum information on <J>(r) one could possibly extract from data. 

Any statistical estimator in general should be able to judge a trial gravitational potential 
<3?(r) and accept it or rule out with a certain significance. Thus, the allowed hypotheses can be 
sorted out at a given confidence level. The virial estimator does not work in this way: it gives 
only a best bet on the mass and a rough guess of its error based on the previous experience with 
numerical simulations of gravitating systems. One of the advantages of the estimator we propose is 
its completeness in the sense that it enables derivation of the errors at any given confidence level. 

The adequate statistical approach becomes especially valuable when one tries to discriminate 
between different models of $(r) that have more than one parameter. Consider for example a set 
of trial models with two parameters p\ and p2- Which model is in best agreement with positions 
and velocities of the observed bodies? What p\ and p2 are allowed at the 90% confidence level? 
The virial relation is not able to answer such questions. It gives instead a relation between p\ 
and p2 (and an approximate error of this relation). The roulette estimator turns out to be much 
more powerful: it constrains both p\ and p2 at any chosen confidence level. We illustrate this 
in § 4 with a simple halo model that has two parameters: the halo size b and mass m. In this 
example, the instantaneous positions and velocities of N test bodies moving in the halo potential 
contain information on m and b which one would like to extract. The roulette method gives a 
best fit (mo, bo) as well as confidence contours on the (m, b) plane, and we test its efficiency with 
Monte-Carlo generated sets of bodies in a known potential (m t ruc> &truc)- We find that the standard 
deviation of (mo, &o) from the true values is as small as 10-15% for N = 32. 

The success of the roulette method is partially due to its ability to use the full 3D information 
on the positions rj and velocities Vj of the test bodies. By contrast, the virial estimator uses 
only absolute values rj and Vi (i = 1,...,N), and thus the angles between Vj and rj are ignored. 
These angles are quite valuable as they contain information on the eccentricities of the orbits. For 
instance, if Vj is parallel to rj it is clear that the body is on a radial (linear) orbit, while the virial 
estimator is not "aware" of that. The advantage of the roulette method in this respect becomes less 
clear in applications where only projected components of rj and Vj are known from observations. 
The method then needs to be modified, which is deferred to a future paper. Already in its 3D 



2 The Bayesian method also makes use of the roulette principle because the satellite distribution function is assumed 
to depend only on the integrals of motion. This form of the distribution function is justified by the strong Jeans 
theorem whose proof relies on the randomness of orbital phases (Appendix 4A in Binney & Tremaine 1987). 
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version, the roulette estimator has useful applications that are briefly discussed in § 5. 



2. Virial estimator 

In this section, we discuss briefly the performance of the virial estimator applied to N test 
bodies orbiting a point mass M. It will be used later as a benchmark. The virial theorem states 
for any bound orbit 

GM 2 

< >=< v z >, (1) 

r 

where angle brackets signify a time average. The same relation holds for a population of N — > oo 
bodies observed at one moment of time if the average is taken over the population. 3 For a set of 
observed positions r; L and velocities V{ one gets the virial mass 

M v converges to M truc at N — > oo, independently of the (unknown) orbital eccentricities of the 
bodies which makes the virial relation tempting for mass estimations. The virial estimator has, 
however, the following drawbacks: 

1. — It does not involve any statistical analysis of the data. Instead of confidence intervals 
the estimator gives a direct formula for M in terms of observed positions and velocities of the test 
bodies. For a finite N, M v is not equal to the true mass, and its mean expectation and variation 
depend on the unknown orbital eccentricities of the bodies e^. 

This is evident in the extreme case of N = 1. Then equation (2) gives M v = G~ 1 v 2 r where v 
and r are the velocity and radial position of the observed body. The mean expectation for M v is 
< M v >= G" 1 < v 2 r > where the average is now taken over the ensemble of random realizations 
of the data set. Each realization can be thought of as a snapshot taken at a random moment of 
time and < v 2 r > equals the time-averaged value of v 2 r for the true orbit (see Appendix), 

< v 2 r >= GM tIUC (l - y) ■ ( 3 ) 

The standard deviation A(v 2 r) can also be calculated, 

A(v 2 r) e /A . 

(4) 



< v 2 r > V2 - e 2 ' 
Thus, both < M v > and A(M„) depend on the unknown e for iV = 1. 



3 The virial relation (1) may be invalid if the population has orbits with strongly varying size a, so that a max /a min 
oo when N — ► oo. 
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A qualitatively similar dependence takes place for N > 1, and therefore it is difficult to derive 
the error of M v without additional assumptions about e^. The error would vanish for circular orbits 
(ej = 0): each body would be expected to give the same M = vfri/G, i = 1, N. In the case of 
large the error can be significant even if N S> 1 and depends on the radial distribution of the 
bodies as we discuss next. 

2. - Test bodies that are at larger radii contribute to the virial estimator with a smaller 
weight and those at the smallest radii make the dominant contribution. If N bodies are spread 
over a significant range of radii, the effective size of the sample is smaller than N, which leads to 
a relatively large statistical error of the estimated mass. 

For illustration, consider a population of orbits with random eccentricity < e < 1 and semi- 
major axis a m i n < a < a max , and make a simple Monte-Carlo simulation. For each body i = 1, N 
we draw randomly e^, Oj, and the orbital time phase of the body. Thus we create a "data set" 
Ti, Vj. We apply equation (2) to the set and obtain M v , which we can compare with the true M. 
Repeating this for many randomly drawn data sets we can study the statistics of M v . 

Figure 1 shows the standard deviation of M v for N = 10, 100, 1000. It depends on the ratio 
a max/o m in : with increasing a m ax/«min the statistical error of M v increases. This happens because 
the bodies at small r (and correspondingly large v) dominate in equation (2), and the data with 
large r is practically lost. 

To avoid this problem one could relate M to < v 2 r > instead of equation (2). This relation, 
however, depends on e^, even at iV — > oo. For example, for bodies on orbits with equal eccentricities 
ei = e one gets (see eq. 3) 

<v 2 r> 

M= G(l-eV2) - (5) 

A projected version of this estimator was proposed by Page (1952) who studied the case of circular 
orbits (e = 0) and by Bahcall & Tremaine (1981) for arbitrary e. The advantage of such estimates 
is that bodies at different r contribute equally, which gives a nice convergence of the estimated M 
with increasing N. Nevertheless, the fact that e is unknown still impedes a precise estimate of M 
even when the exact < v 2 r > is known (N — > oo). 

3. — The virial estimator uses only absolute values r« and v j and ignores the angles between 
Vj and Ti. 

Finally, we note that the virial estimator is inefficient when applied to non-Keplerian potentials. 
The virial relation for a general potential is given by 

N 

2(v? + ri .fi)=0, (6) 

i=i 



where fj are the gravitational accelerations of the test bodies. It gives one condition and, if the 
potential has more than one unknown parameter, the relation is not able to determine them. 
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Fig. 1. — Standard deviation of the virial estimate M v for a point mass with N observed satellites. 
The orbits of test bodies are assumed to have random orbital eccentricities < e < 1 and semi- 
major axes a m j n < a < a max . 



- 7- 



3. Orbital roulette 

Suppose we have measured a current 3D position r and a velocity v of a satellite orbiting 
a central massive object. What can we say about the central mass M based on this "snapshot" 
information only? First of all, assuming that the satellite is on a bound orbit (its age is longer than 
the orbital period), we get a lower bound M m - m = v 2 r/2G. Then we note that if M is very large, 
GM/r 3> v 2 /2, then the satellite has to be extremely close to its apocenter, and if M is small, 
M — > M m i n , it has to be near the pericenter (Fig. 2). The probability of either extreme is low 
because the snapshot is taken at a random moment of time and, in most cases, the satellite should 
be somewhere between the pericenter and apocenter. Thus, the requirement of a random orbital 
time phase can constrain M from above and below at a given confidence level. If we have a snapshot 
of iV 3> 1 satellites, the constraint becomes tight. Below we quantify this constraint and develop the 
mass estimator based on the random-phase principle. The simplest way to check the randomness 
is by comparing the mean phase with its expected value (§ 3.1) and a more sophisticated method 
checks the whole distribution of N phases (§ 3.2). 



3.1. Mean orbital phase 

Let us denote the unknown orbital period of a satellite by T and the time since last passage of 
the pericenter by t p . Since we know the current r and v, the whole orbit and t p are unambiguously 
calculated for any given M tr i a i- For a known orbit, we need a measure of how close the satellite is 
to the pericenter. We define this measure as time separating the satellite from the nearest passage 
of the pericenter (in the past or the future) and normalize it by T/2, 

so that g can take values between and 1. This definition takes into account that the pericenter 
is closer in the past if the current radial component of velocity v r > 0, and in the future if v r < 0. 
The limit M tria i — ► M m ; n gives g — > (pericenter), and M tria i — > oo gives g — > 1 (apocenter). For 
the true M and a random snapshot time, the expected g obeys Poisson statistics: it has a flat 
probability distribution between and 1 with the mean expectation value < g >= 1/2 and the 
standard deviation Ag = 12 -1 / 2 . 

Now suppose we have a snapshot of N 3> 1 satellites with measured rj and Vj, i = 1, N. For 
a given M tria i we can calculate the orbit of each satellite and the corresponding 5i(M tria i). Then 
we define the mean phase, 

1 N 

S(Mtrial) = & ( M trial)- (8) 

i=l 

Again, g — > for small M tr i a i and g — > 1 for large M tr i a i- By contrast, for the true M, the expected 
g is described by a narrow probability distribution f(g) which peaks at 1/2. At N S> 1, this 
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Fig. 2. — Orbit reconstruction for a body with a given position and velocity vector v. The recon- 
structed orbit depends on the assumed central mass M. A small M -► M min = rv 2 /2G gives a 
large orbit and then the current position is much closer to the pericenter (shown by the open circle) 
than a distant apocenter far outside the figure. A large M (100M m i n in this example) gives a small 
orbit and places the body almost exactly at the apocenter. 
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distribution is Gaussian according to the central limit theorem, 



/( ^ ) = 71^ exp 



2a 2 

So, one can find the gravitating mass by adjusting M so that 



a = ^=. (9) 
Vl2N 



g(M) = i±(12iV)- 1 /2. (10 ) 

The estimate is well defined since g is a monotonic function of M. The width a of the distribution (9) 
sets the error of the estimate. Expanding equation (10) near the best-fit point g = 0.5, one gets 
the error in M, 

dg 



™ = (12iV)-V 2 
M V ' 



M- 



dM 



(11) 

9=1/2 



M(dg/dM) = O(l) for any large N (variation of M by a factor ~ 2 around the best-fit value 
induces a change of g comparable to 1/2) and hence AM/M = 0{N^ 1 ^ 2 ). The estimate of M 
defined by equation (10) converges to M true at N — ► oo. 

The analysis of orbital phases gi is similar to testing a roulette wheel. For a fair roulette, 
the ball angular position on the wheel must be random. Suppose we have a reference position, 
e.g. Zero, and make N experiments with the roulette. Then we get N random angular deviations 
from Zero, < ai < 180°, i = 1, N. If the wheel is unbiased, one expects ai to follow Poisson 
statistics with a mean value di = [0.5 ± (12A0~ 1/2 ] x 180°. There is a full analogy with our problem 
as the magnitude g = a/180° obeys the same probability distribution (eq. 9). The only difference is 
that the random variable of orbital roulette is the time phase of a test body rather than its angular 
position. 

Equation (10) gives a best-bet value of M and an estimate of its error. We now aim to obtain 
a more complete solution to the problem: the allowed interval for M at a given confidence level C. 
Mathematically, this is formulated as follows. For a given number < ( < 1/2 we evaluate M + such 
that the probability of M true > M + equals £, and M_ such that the probability of M true < M_ 
equals £. The interval M_(£) < M < M + (£) is the mass measurement at the confidence level 
C = 1 - 2£. 

Let us define cumulative distributions 

P-(g)= I" f(g)dg, P+(g)= f fm~g- (12) 

JO J g 

P- + P + = 1 for any M tria i. Then define <?-(£) and <?+(£) so that 

P-{g-)=i, P+{g+)=i. (13) 

Note that g± are unique functions of £ (they depend only on N as a parameter) and the snapshot 
data i-j and Vj did not appear in the definitions (12,13). The data determine g(M) and appear in 
the final equation for M±(£), 

g(M-)=g-, g(M + )=g + . (14) 
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One can prove that M± defined in this way give the correct confidence interval for M at the 
confidence level C = 1 — 2£. Consider the ensemble of all possible snapshots of our test bodies 
taken at random moments of time. For each snapshot one has M_(£) defined in equation (14). 
What is the fraction of "bad" snapshots with M truc < M_(£), outside the confidence interval? For 
these snapshots, 

fftrue = g(M true ) <g-{0, (15) 

[this is equivalent to M truc < M_(£) because g(M) is a monotonic function for any snapshot]. 
The probability to get a snapshot with g true < g_ is P-(g-), and it equals £ by definition of g_. 
Hence, the probability of a snapshot with M tTUC < M_(£) equals £. Analagously, one shows that 
the probability of a snapshot with M truc > M+(£) equals £. 

P-(g) ^ 1 signals M tria i < M truc , and P+{g) <C 1 signals M tria i > M true . The best bet for 
M is defined by P_(M) = P+(M) = 1/2, which is equivalent to g(M) = 0.5. The best bet is 
supplemented with the confidence intervals (M_,M + ) at any confidence level, which also gives an 
explicit probability distribution of M, 

*< M > = m = - w = H /te[M »- (16> 

This probability distribution is the measurement of M with the mean-phase method. 



3.2. Cumulative distribution of orbital phases 

When the trial mass is close to M true , the inferred phases gi should be consistent with a Poisson 
distribution, which can be used as a test for M tr i a i. Such a test has the advantage of using the 
whole distribution of g± rather than just the average g. Consider, for example, gi half of which 
equal and the other half equal 1. This is clearly inconsistent with a Poisson distribution, yet g 
has the right value of 1/2 and the mean-phase method will not notice the inconsistency. 

Testing a distribution for consistency with an expected/model distribution (null hypothesis) is a 
standard problem of mathematical statistics. It is done by comparing the corresponding cumulative 
distributions. (Kolmogorov 1941). The cumulative distribution of a given set g^ (i = 1,...,N) is a 
stepped function, < F(g) < 1, that increases by 1/N at each <?j. The null hypothesis in our case 
is a Poisson process, which gives the mean expectation < F(g) >= g at any g (and for any N). 
Random realizations of the Poisson set gi have F{g) fluctuating around < F(g) > with a binomial 
distribution, 

= (FN)KN- FNV 9FN{1 - S) "™ <17> 
<F(g)>=g, VarfFfe)] = A'[F(g)] = (18) 

where Var[X] and A[X] denote variance and standard deviation of X. The unknown mass M should 
be adjusted so that the cumulative distribution of gi(M) is consistent with a Poisson process. 
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One then needs a measure of the difference between the obtained F(g) (which depends on M 
as a parameter) and the mean expectation < F(g) >= g. The simplest measure would be the 
Kolmogorov-Smirnov test: one defines Vks = max \F(g) — g\, where the maximum is searched over 
the whole interval < g < 1 (Kolmogorov 1941). The probability distribution of Vks 1S known for 
a true Poisson process (e.g., Press et al. 1992) and getting the actual Vks f ar i n the tail of this 
distribution signals a low probability of consistency with the null hypothesis. 

The Kolmogorov-Smirnov measure is most sensitive to the values of F(g) at median g ~ 0.5 
where Yai[F(g)} is largest (see eq. 18). In our case, a measure sensitive at the tails g « and g pt 1 
would be preferable because, when M tria i deviates from M true , gi cluster near or 1. We therefore 
use a different measure 

v= f 1 [f (a)- < f (g) >1 2 d jy /' l F (s) - s? d (19) 

J V«r[F(s)] J, j(l-j) V ' 

It gives more weight to the tails and is known as the Anderson-Darling test (Anderson & Darling 
1952). It can be viewed as x 2 f° r a fit of F(g) by F = g. Equation (19) becomes the usual formula 
for x 2 if the integral is replaced by a finite sum over dg = 1/K with K — > oo. 

The probability distribution of V for a true Poisson process can be calculated numerically by 
the Monte-Carlo technique. We randomly generate many Poisson sets gi [i = 1,...,N), calculate 
V for each set, and find the histogram of V. The found probability distribution p(V) is shown in 
Figure 3. It does not depend on N as long as N S> 1. At V > 2 (p < 0.1) the distribution is very 
well approximated by a simple empirical formula, 

p(V) = W- v/2 , V>2. (20) 

Having gi(M tria ]) and their V, one can quantify the inconsistency of the trial with a Poisson 
process. The significance level of inconsistency, 

[•OO 

= / P (V')dV', (21) 

is the probability that V computed for a random Poisson set g[ is greater than V. A low £ means 
that it is unlikely to draw occasionally a Poisson set g[ with V' as large as our V, and so our M tria i 
is rejected at the significance level £. 

For V(M) > 2, one can use equation (20) to get a simple explicit formula 

C{V) = mT0 10 ~ y/2 ' V>2 ' (22) 
All M tr iaj that give V(M tr - ia x) > — 2 log 10 [(£/2) In 10]) are rejected at the significance level £. 

The function T^(Mt r i a i) may be non-monotonic (especially near its minimum), and a number 
k > 2 of mass intervals might be rejected. M true belongs to the remaining regions at the confidence 
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level C = 1 — k£. In practice, however, our Monte-Carlo simulations show that the case k > 2 
rarely occurs for £ of practical interest, and we count only two rejected intervals that cover the 
tails of low mass (M m i n , M_) and high mass (M+, oo). M true belongs to the interval (M_, M + ) at 
confidence level C = 1 — 2£. 



3.3. Practical algorithm 

We summarize now the practical steps to derive the central point mass M from data. Given 
a measured position r and velocity v of a test body, one first determines its orbital phase g as a 
function of M. An explicit formula for g is derived in Appendix A, 

= -(2^-esin2V>), (23) 

7T 

where 

and the orbital energy E and eccentricity e are expressed in terms of r, v, and unknown M in 
Appendix A (eqs. 33,34). Equation (23) gives gi(M) for each test body % = 1,...,N. 

The mean-phase variant of the roulette estimator is simplest to apply. It defines the best-fit M 
by equation g(M) = 9i/N = 1/2. Since g increases monotonically with M, only one solution to 
equation g(M) = 1/2 can exist. The solution is easily found with the bisection method. The initial 
interval (M m i n ,M max ) for bisection can be chosen as follows. M m \ D is given by the condition that 
all N orbits are bound: M m \ n = max.(yf /2Gri). M max should be taken sufficiently large, so that 
g(M max ) > 1/2. Note that no solution exists if g(M m [ n ) > 1/2. It may happen in rare cases where 
the observed bodies accidentally cluster near their true apocenters and g(M truc ) is large. Such 
snapshots may have g(M m i n ) > 1/2, and then the best-fit mass is not defined by the mean-phase 
estimator. Only confidence intervals should be considered in this case. 

The confidence intervals provide the most useful information. They show the allowed interval 
for M at a given confidence level C = 1 — 2£ where £ can be chosen at any small value. In 
practice, measurements at the 90% confidence level are often considered, which corresponds to 
£ = 0.05. Confidence intervals (M_,M + ) derived with the mean-phase estimator are given by 
equation (14) where <7±(£) are defined by equations (12,13). Note that f(g) is precisely Gaussian 
only for N S> 1, and at small N one should use a different formula: the convolution of N constant 
functions fi(g) = 1, < g < 1 (the convolution approaches Gaussian quickly with increasing N; 
e.g. at N = 5 the Gaussian approximation is already good). 

The test for the cumulative distribution of orbital phases is slightly more complicated, however, 
it has the advantage of utilizing all the available data. For any trial mass M one can calculate the 
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Fig. 3. — Distribution of V (eq. 19) for a true Poisson process with N » 1. The solid curve 
shows the result of numerical (Monte-Carlo) calculation. The dashed line shows the approximation 
p(V) = l(T y / 2 , which is valid at V > 2 (p < 0.1). 
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Anderson-Darling V(M) (eq. 19), which can be written as 



V 



N 
i=0 



N 



9% 



1 



- In 

N 



1-9, 



»+i 



9i 



+ 9i- 9i+i 



(25) 



Formal go = and g^+i = 1 ar e used in this expression. The best-fit M is such that V(M) is 
minimal. This minimum can be found numerically. Alternatively, one solves the algebraic equation 
dV/dM = 0, which can be written down analytically. The derivative dV/dM is 

N H " r 2 (iv-i) + i (2,-i)l (26) 



dV 



dM 



1 dgi 



N 2 ^ dM 

i=l 



9i 



9i 



and the expression for dg/dM is derived in Appendix A. This gives an explicit analytical form to 
equation ^{M) = 0. The equation may, however, have many roots, i.e., V(M) may have many 
local minima and maxima around the global minimum. Therefore in practice it is easier to find the 
minimum numerically. 

Once V(M) is evaluated, one should find the intervals of M that can be rejected at a given 
significance level £ as explained in the end of § 3.2. If £ < 0.1 (i.e. the estimate is done on a 
confidence level C = 1 — 2£ > 80%), the analytical approximation to V(£) can be used. Then the 
boundaries of the rejected intervals are solutions of the equation V(M) = — 21og 10 [(£/2) In 10]. 



3.4. Handling data uncertainties 

The described method is easy to apply to a snapshot of N bodies with precisely measured rj 
and Vj. Real data, however, have uncertainties. Suppose a measurement gives 3D distributions 
p(vj) and p{ri) rather than numbers Vj and r^; for example, they may be Gaussian distributions 
whose widths represent the measurement errors of rj and Vj. Then, for a given M tr i a i, the orbital 
parameters of the bodies are also described by probability distributions. As a result one obtains N 
probability distributions Pi(gi) and p(g) instead of numbers gi and g. 

The p(g) found for M tr ; a i should be compared with the expected distribution f(g) (eq. 9) 
and their consistency should be evaluated. This can be done using the Kolmogorov-Smirnov or 
Anderson-Darling tests. In practice, it is easier to use the Monte-Carlo technique that incorporates 
the data uncertainties directly in the roulette test for M tr i a i- In each Monte-Carlo event, one draws 
randomly Vj and from the measured distributions p(vj) and p{v i ), calculates the orbits, and finds 
the mean orbital phase g of the N satellites. Then one compares g with g' for N randomly drawn 
numbers < g[ < 1. Repeating this comparison for many Monte-Carlo Vj, rj and g' , one finds 
the probability P_(M tria i) that g > g' and the corresponding P + (M tria i) = 1 — P_. Thus, one gets 
the same estimator as the one described in § 3.1 but now it allows for the data uncertainties. The 
best-fit mass is found from the condition P_ = P + = 1/2 and the confidence intervals are found 
from P_ = P + = £. In the case of large data uncertainties, the functions P±(M tria i) are modified 
and the range of M consistent with the data becomes larger. 
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The mass estimator based on the Anderson-Darling test of the cumulative distribution of gi is 
modified in a similar way to incorporate the data uncertainties. 

In a separate paper on the Galactic Center, we present an analysis of a real data set and 
demonstrate how the measurement errors are taken into account. 

3.5. Numerical tests 

We now check the performance of the roulette estimator by direct Monte-Carlo simulations. 
With this purpose, we have developed a numerical code that calculates the best-fit mass and 
confidence intervals with both mean-phase and cumulative-distribution variants of the estimator. 

First of all we check that the confidence intervals are correct. We take N bodies orbiting a 
central mass M true with randomly chosen eccentricities and semi-major axes. Each body is taken 
at a random moment of time. This gives us a simulated data set to which we apply the estimator 
and obtain the 90% confidence intervals (M_,M + ). We repeat this procedure for many generated 
data sets and count the number of cases where M true is outside our confidence interval. These cases 
should make 10% and we have checked that this is indeed so. 

Then we compare the roulette method with the virial estimator. Since the virial estimator does 
not give confidence intervals, we have to restrict the comparison to the best-fit values of M. We 
denote them by M r and M v for the roulette and virial estimators, respectively. Computer-generated 
random data sets allow us to study the statistics of M r and M v . 

The advantage of the roulette estimator becomes significant when the generated orbits have 
significantly different sizes. For illustration, we draw random semi-major axes of the orbits from 
an interval < a < a max . We also draw a random eccentricity < e < 1 for each orbit. Each 
generated data set of N bodies is analyzed with three methods: (1) virial, (2) roulette mean-phase 
test, and (3) roulette cumulative-distribution test. The results of analysis of 10 5 data sets are 
shown in Figure 4. We observe that for N > 4 the roulette method gives a more precise estimate of 
M than the virial estimator and shows a much faster convergence to M truc with increasing N. We 
also observe that the mean-phase estimator works practically as well as the more complete analysis 
of the cumulative distribution. 

The roulette method is illustrated in Figure 5 for six randomly selected data sets with N = 10. 
The best-fits and 90% confidence intervals obtained with the two roulette estimators strongly 
correlate with each other. 
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Fig. 4. — Performance of the virial and roulette estimators tested with N bodies on orbits with 
random eccentricities < e < 1 and semi-major axes < a < a max - The mean value and standard 
deviation of the estimated mass are shown in the figure as a function of N. Two versions of the 
roulette estimator are shown: the mean-phase method and the Anderson-Darling (AD) test for the 
cumulative distribution of phases. 
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Fig. 5. — Roulette analysis of six randomly selected data sets with N = 10 from the sample in 
Figure 4. Each panel displays the results for one data set. The measure V of the cumulative- 
distribution estimator is found as a function of trial mass M and shown by solid curve. The 
minimum of V(M) gives the best-fit M, and V = 2.48 defines the 90% confidence interval for 
M (shown by the horizontal line). The mean-phase estimator compares g{M) with 1/2, and the 
dashed curve shows h = |5-0.5|(12iV) 1 / 2 as a function of trial M. Then the best-fit is where h = 0, 
and the 90% confidence interval is defined by h = 1.643 (0.35 < g < 0.65). 
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4. The case of distributed mass 

Suppose we observe N bodies moving in an unknown potential <J>(r) which is created by a 
distributed mass with an unknown density profile. Hereafter for brevity we call the gravitating 
mass "halo" (keeping in mind a dark-matter halo of a galaxy) but it could also be a stellar cluster 
or a cluster of galaxies. We need only assume that the cluster has a sufficient number of members, 
so that gravity fluctuations due to motions in the cluster average out, a quasi-steady <l?(r) is well 
defined, and each observed body moves on a well-defined bound orbit in this potential. We now 
apply the roulette method to estimate <E>(r). 

We limit our consideration to spherically symmetric potentials; then an orbital motion is always 
confined to a plane. Bound orbits in a potential 3>(r) are not closed in general, however, they still 
have well-defined pericenter n and apocenter r^. The bodies move periodically from the pericenter 
to the apocenter and back with a period T. This "radial" period is in general not equal to the 
period of the azimuthal motion. A typical example of such an orbit is a rosette: a superposition of 
radial oscillations and azimuthal precession (e.g. Binney & Tremaine, 1987). 

Since r\ and r2 are well defined for each orbit, a moving body has a well defined orbital phase 
g: the time of motion from the current radius r to the nearest pericenter (in the past or in the 
future) divided by half period, T/2. If the body is observed at a random moment of time then the 
phase should be a random number < g < 1 and we can apply the roulette test. 

Any trial model $(r) can be tested as it gives certain phases gi to the observed bodies i = 
1, N. The model is good if gi are consistent with the Poisson distribution, which can be checked 
in exactly the same way as it was done for point-mass potentials in § 3. In this way, the roulette 
test sorts out good models <3?(r) from bad ones. 

In practice, one deals with a parametrized family of potentials assuming a certain functional 
shape of $>(r). Estimating <3>(r) is then reduced to constraining the allowed range of a small 
number of parameters. The parametrization is consistent with the data if it passes the roulette 
test with some set of parameters. The test will signal if the assumed parametrization of 3>(r) is far 
from reality: the model will not be able to produce Poisson phases if N is sufficiently large. For 
example, a point-mass model with one parameter M will not give Poisson g, L if the bodies actually 
move inside a homogeneous halo. Thus, the roulette method checks the assumed parametrization 
and simultaneously finds the best-fit parameters. In concrete problems, one may have an idea of 
the possible functional shape of the potential. For example, dark-matter halos are expected to have 
certain shapes predicted by cosmological simulations of structure formation (e.g. Navarro, Frenk, 
k White 1996). 

For the illustrative purpose, we consider below a simple potential with two parameters b and 

m, 

*w=- 6+( ,f;v - < 27 > 
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This is the well-known isochron potential (see e.g. Binney & Tremaine 1987); it describes a gravi- 
tating mass m distributed in space so that most of it resides inside the characteristic radius b. The 
corresponding mass density is given by the Poisson equation div <I> = Air p. At r <C b, the density is 
p const = (3/167r)(m/& 3 ). At r > b the density falls off quickly and <J?(r) tends to the point-mass 
potential <J?(r) = —Gm/r. The equation of motion in the isochron potential can be integrated 
analytically and the orbit shape can be found (see Appendix B). 

4.1. Phase-energy correlation 

Let us test the efficiency of the roulette estimator with Monte-Carlo generated sets of N test 
bodies moving in an isochron potential. The true potential has known parameters m truc and &true- 
We assume in this test that the bodies have random orbital energy Eq < E < and random 
eccentricity < e < e max defined in Appendix. 4 The minimum Eq = — (Gm true /26 tr ue) = ^true(O) 
corresponds to a body at rest at r = 0, and e max (E) = 1 — E/Eq is the maximum eccentricity for 
an orbit with a given E. 

The created data set is studied by the roulette estimator which does not know mtrue, btrue- 

The 

estimator calculates the orbital phases gi(m,b) and compares them with the Poisson distribution 
using the Anderson-Darling measure V as described in § 3. Then it finds the best fit (mo, &o) where 
V(m, b) has a minimum. 

We, however, know m tr ue and 6 true and can check the accuracy of the obtained estimate (mo, bo). 
This check is made in Figure 6 (upper panel) for 300 simulated data sets with N = 32. We plot 
the results in the q-m plane, where q = m}/ 3 /b is a density parameter. The choice of q and m as 
two independent parameters (instead of b and m) has a reason that will become clear below. 

One can see relatively large deviations of (qo,mo) from (<?true, WHrue)- Interestingly, the devia- 
tions are concentrated in an elongated region in the q — m plane and their origin can be understood. 
First, we note that a wrong choice of m affects mostly gi of bodies at r > b, and a wrong choice of 
q affects gi of bodies at r < b. Consider, for example, m > m truc . The orbital phases of the bodies 
at r > b are then biased to the apocenter. This can be compensated if we choose q < (?true : then 
the bodies at r < b will be biased to the pericenter. Thus, m > m truc and q < (/true can give the 
correct g = 1/2 even when they deviate significantly from m^ue and (ftrue- The total distribution 
of gi is then biased to both apocenter and pericenter, and there is a lack of intermediate gi ~ 1/2. 
At large N, the Anderson-Darling test notices such a non-Poisson distribution of gi and rules out 
(q,m) that are far from (<ftrue, ^true)- With increasing N, the acceptable (q,m) converge to the 
true values. However, this convergence is relatively slow. 

This problem is general for estimations of distributed-mass potentials that have more than 



4 More precisely, we assume in this example a Poisson probability distribution for \J— E between and \J—Eo and 
a Poisson distribution for e between and e max . 
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Fig. 6. — 300 best fits to the isochron potential probed with N = 32 bodies randomly drawn from 
a population with random orbital energy and eccentricity (see the text). The fits are shown in the 
q — m plane where q = m^'^/b; m and b are the mass and size of the isochron halo (eq. 27). (a) 
Roulette method (AD test), (b) Casino method. 10% worst points are in red, 50% best points are 
in blue, and 40% intermediate points are in green. 



- 21 - 



one parameter: there are directions in the parameter space where g = 1/2 is satisfied and the 
non-flatness of the g± distribution has to be detected in order to rule out false parameters. The 
non-flatness is more difficult to detect than a deviation of g from 1/2 because of a lower signal-to- 
noise ratio, and therefore larger N are required to estimate the potential accurately. 

Fortunately, there is an efficient way to avoid this problem and improve the accuracy of the 
roulette estimator. Indeed, we note that the large deviations of (<7o> m o) from ((ftrue> mt ruc ) are 
accompanied by a strong intrinsic correlation between E; L and gi. In our example above, wiq > m trU e 
and qo < (/true, the bodies on small orbits (small E) are found near the pericenter and the bodies 
on large orbits (large E) are near the apocenter. The orbital energy Ei is an unknown integral 
of motion (the other integral, angular momentum lj = x Vj, is known from the data), and our 
observation can be generalized as follows: large deviations of the best-fit parameters from the true 
values occur when the orbital phases and unknown integrals of motion are adjusted in a special way 
and show an intrinsic correlation. Real, truly random, gi should show no correlation with integrals 
of motion. Thus, the roulette method can be improved by tracking this correlation and making 
sure it is absent for the best fits. In particular, in our halo problem, we should track the E — g 
correlation. 

A standard way to detect a monotonic (positive or negative) correlation is to calculate the 
cross-correlation coefficient XX ~~ E)(9i ~ <?)• We will use a better measure that would signal any 
type of correlation between gi and Ei . First we express the problem mathematically in a convenient 
way. Given gi and Ei we can sort the bodies in two ways: by g and by E. Let us denote their 
phase and energy numbers by i and j, so that gi + \ > gi and Ej+i > Ej. The absence of correlation 
between E and g means that the permutation z — > j is random. This permutation can be viewed 
as a mixing of ./V cards (bodies), and in a fair card game we expect a truly random mixing. Plotted 
on the i — j plane, a random permutation is represented by a set of ./V points that should be 
consistent with a homogeneous 2D distribution in the square (1, N) x (1, N). Inhomogeneity would 
signal a correlation. 

The consistency with homogeneity can be checked as follows. Consider one point It 
divides the square into four quadrants: (1) i' < i, j' < j, (2) i' < i, j' > j, (3) i' > i, j' < j, (4) 
i' > i, j' > j. Denote the number of points that fall into each quadrant as Ni, N2, N3, N4. For a 
random mixing the mean expectation values for these numbers are proportional to the areas of 
the corresponding quadrants. Only one of the four numbers is independent, which we choose to be 
Ni (we have N 2 = i - 1 - N 1: N 3 = j - 1 - Aq, and N 4 = N-i-N 3 = N-i-j + l + N{). Random 
mixings give the hypergeometric probability distribution for N\ with the mean expectation 
value and variance, 



< Aq >= (j - 1) 



(i-1) 

(N-iy 



Var(Aq) 



(N-i) 
(N-2) 



(j-l)(N-j) 
(N - l) 2 



(28) 



Aq is consistent with random mixing if its deviation from < Aq > is comparable to [Var(Aq)] 1 / 2 . 
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The inconsistency is measured by 

" = %^) >>2 - < 29 > 

Finally we define a measure that should signal deviations from a homogeneous distribution in any 
part of the square (l,N) x (l,N), 

N-l 

V c =^V c [i,j(i)]. (30) 

V c plays the same role for the card test [randomness of mixing as the Anderson-Darling 
measure V played for the roulette test (randomness of phase Using the Monte-Carlo technique 
we have calculated numerically the statistics of V c for truly random card mixings. We thus have 
the probability distribution of V c , its mean expectation < V c >, and variance \ar(V c ). A mixing 

with some V c will not pass our card test at a significance level £ if the probability that a fair 
V' c > V c equals £. 



4.2. Casino 

For the true potential $trueM we should have random orbital phases <j>j and random energy- 
phase mixing Thus we expect to deal with a fair casino, 

fair casino = fair roulette + fair cards. (31) 

For a false trial potentials <J>(r) we will detect that its casino is unfair by looking at the combination, 

^V^+V^W <32> 

Our best bet for <J?(r) is the potential that minimizes x 2 ■ 

A fair casino has a certain cumulative probability distribution P(x 2 ) which we have calculated 
numerically. It is shown in Figure 7. For example, one can see from the figure that a model that gives 
X 2 ~ 16 is good with 90% confidence (10% significance of inconsistency with the null hypothesis of 
fair casino). Using P(x 2 ) one can sort out acceptable potentials at any given confidence level. 

Let us illustrate the efficiency of the casino method by applying it to the isochron potential. 
We repeat our Monte-Carlo simulation shown in the upper panel of Figure 6, but now we get the 
best fit by minimizing x 2 rather than V. The results are shown in the lower panel of Figure 6. 
The accuracy of the casino estimator is improved significantly compared with the simple roulette. 
The standard deviations of the best fits from (</true> ^true) are now 11-13%, and 90% of the best 
fits show deviations less than ~ 25%. 



We think that the casino method extracts completely the available information from the data 
and does the best estimate of <3? that could possibly be done. As an illustration, we show in Figure 8 
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Fig. 7. — Cumulative distribution of the casino measure x 2 defined in equation (32). Short-dashed, 
long-dashed, and solid curves show the cases N = fO, 32, and fOO, respectively. 



-24- 



the results of the method application to three Monte-Carlo-generated data sets with N=100. In a 
similar way, it would be applied to real data. The only difference is that here we know the true 
potential (and the method does not). The results are presented as 90% confidence area on the q — m 
plane. The data set analyzed in panel (a) is drawn from a population with random orbital energy 
E between and E = — (Gmt rU e/26 true ) and random e. It gives a compact confidence area and 
both q and m are correctly measured with accuracy better than 10%. Panel (b) shows the result for 
a data set with E randomly drawn from the interval Eq < E < (19/20) .Eo- In this case, the bodies 
sit well inside b and the method cannot say anything about the total mass of the halo m. This is 
evident for us because we know mt rue and &tme> so we know that the bodies are inside b. However, 
the method initially did not know that. It finds that it cannot say anything about m and hence 
finds that the observed bodies move in a central part of a distributed mass. It measures well (with 
4% accuracy) all it could possibly measure: the central density p = (3/167r)g 3 . Panel (c) shows the 
case where E is drawn from the interval Eq/20 < E < 0. Here the dominant majority of the bodies 
are far outside b and the method measures well the total mass m but cannot say anything about b 
and p. Thus it finds that the test bodies move in essentially a point-mass potential and estimates 
this mass with accuracy better than 10%. 

We have also checked the method ability to reject incorrect parametrizations of 3>(r). For 
example, if the test bodies move inside a homogeneous halo, N = 10 — 15 is enough to reject 
the point-mass model: the best-fit then has too high \ 2 an d does not pass the casino test. The 
confidence level of rejection increases quickly with N. 

5. Discussion 
5.1. Limitations of the roulette method 

There are two situations where the roulette method should be applied with caution. 

1. — - Orbital phases of the test bodies are not independent. This happens if, for example, the 
satellites are in orbital resonances with each other and the phase-locking effect takes place, as is 
the case for planets in the outer solar system. 

2. — - The snapshot of the test bodies has a finite size comparable to the size of their orbits. Then 
there is a risk that the orbits extend beyond the snapshot boundary, and some bodies are observed 
preferentially near their pericenters just because they would not be seen otherwise. This bias to the 
pericenter can be significant only for satellites with highly eccentric orbits and it can be corrected: 
for each orbit reconstructed with a trial potential, one can check whether it extends beyond the 
snapshot boundary and calculate the fraction of the orbital period that the body spends beyond 
the boundary. 
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5.2. Applications 

We briefly discuss below three astrophysical situations where the proposed method can be 
applied in its full 3D version developed in this paper. 5 

5.2.1. The Oort limit 

The mass content of the Galactic disk in the solar neighborhood is inferred from the vertical 
motions of the nearby stars with respect to the galactic plane (see, e.g., chapt. 4.2 in Binney & 
Tremaine 1987). Mathematically, the problem reduces to a reconstruction of a one-dimensional 
potential <&(z) from a snapshot of one-dimensional motions of N test bodies. 

The Jeans equation was previously used to reconstruct $>(z), which is not the most efficient 
method because it involves (arbitrary) binning of the data and numerical triple differentiation. 
Orbital roulette, in its full casino version, should be able to use the data optimally. It is expected 
to find &(z) that gives the stars random phases of vertical motion (roulette test), and the phases 
should be uncorrelated with the stars' orbital energies (card test). These conditions are necessary 
and sufficient for &(z) to be realistic. 

5.2.2. The mass of Sgr A* 

Observations of stellar motions in the Galactic Center allow one to estimate the mass of the 
central dark object Sgr A*, which is believed to be a giant black hole (Genzel et. al. 2000; Schodel 
et. al. 2002, Ghez et. al. 2003). Normally, only sky-projected positions of stars are available, with 
the exception of a few tightly bound stars whose orbits have been partially mapped out. 

By analyzing the data of Genzel et. al. (2000) we have recently found that the young stellar 
population in the Galactic center forms a thin disk (Levin & Beloborodov 2003). The inferred 
orientation of the disk plane was confirmed by further analysis of Genzel et. al. (2003). This 
finding gives a full 3D information for instantaneous motions of ten disk stars, which can be used 
to estimate the gravitational potential. We have applied the orbital roulette to these data and 
obtained a new independent estimate of the black hole mass (in preparation). 

5.2.3. Galactic potential probed by SIM 

Observations of satellites of our Galaxy constrain its gravitational potential and mass content. 
Within a decade, the Space Interferometry Mission (SIM) can provide the full 3-D information on 



5 We are grateful to Scott Tremaine who pointed out the applications outlined in §§ 5.2.1 and 5.2.3. 
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the motion of tens of the satellites, and the data should be used most efficiently. 

The most advanced method to date, which is based on the Bayesian analysis, was developed 
by Little & Tremaine (1987) and Kochanek (1996). Little & Tremaine (1987) assumed that the 
satellites' velocities are either radial or have isotropic distribution while Kochanek (1996) assumed 
some a priori form of the satellites' distribution function. The orbital roulette, combined with 
the 3-D data, will alleviate the need for a priori assumptions. In the previous section we have 
demonstrated how our method estimates the mass and size of a spherically symmetric halo. A 
similar analysis can be done for the SIM data with a realistic model of the Galactic potential. 

We thank Scott Tremaine for many helpful discussions. We also thank Mark Wilkinson, 
the referee, for his suggestions that helped improve the presentation of the paper. Both authors 
acknowledge financial support from NSERC. A.M.B. was supported by Alfred P. Sloan Foundation. 



Appendix A: Keplerian orbits 

A Keplerian orbit around point mass M is characterized by two integrals of motion: orbital 
energy and angular momentum, 



v 



2 



GM 



E = , 1 = r x v. (33) 

2 r 

The semi-major axis a and eccentricity e of the orbit are given by 

GM 9 2\E\l 2 , N 

a= Wy i-'-mn- <34> 

The pericenter and apocenter radii are 

n=a(l — e), r 2 = a(l + e). (35) 

For the orbital roulette we need to determine the time of motion from n to a given r. It is 
convenient to calculate the time using the relation 



dr 



1/2 



dt = ~ = '\W\i rd ^' (36) 



where ip is defined by 



- 1 * = J ^ <-) 

so that r = r\ at i/j = and r = r 2 at ip = ir/2. The full orbital period corresponds to ijj changing 
from to it. The time of motion from a given r to T\ is 

, , T ( , sin2^\ . , 

t v { T ) = -U-e—f-\, (38) 
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where 



T = 2t{ n ^r 2 ) 



itGM 



V2\E\ 3 / 2 ' 

is the orbital period. Time average of a magnitude Z over the orbit is given by 



< Z >-- 



-rj: 



Zdt=-[ — 
T \ E 



1/2 r n/2 



Zrdip. 



(39) 



(40) 



Suppose an observed body has a measured position r and velocity v, and denote the angle 
between r and v by a. If we do not know the central mass M, the orbital parameters of the body 
and its phase g = t p (r)(T/2)~ 1 will depend on the assumed M or x = GM/v 2 r. In particular, we 
have 

(2x - 1) 



e(M) 



1 



x^ 



■ sin 2 a 



1 — xe 



1/2 



2xe 



dg 



dM ireMx 



sin 2 ip(M) 

ff (M) = -^-|sin2^), 
(1 — e cos 2iji) 



sin 2ip 



\ 2 9 

1 \ sin a 



x 



sin a sin 2ip 



(41) 

(42) 
(43) 

(44) 



These expressions are used in section 3.5. 



Appendix B: Isochron orbits 



In contrast to the Keplerian case, we now have a characteristic scale b. It is convenient to use 
a dimensionless variable s instead of r, 



s = l + yi + p, r = by / s{s-2). 
The pericenter and apocenter radii are derived from the equation, 

..2 p 



JL = E-$ = 0, 

2 2r 2 



which gives 



1 , Eo 
ai>2 = l + -=F 



1/2 



(45) 



(46) 



(47) 



Here -Eo = ^(O) = — (Gm/2b) is the minimum possible energy for an orbit in the isochron potential. 
It is convenient to define a "semi-major" axis a in the s-space as 

Gm 



2a = b(si + s 2 -2) 



\E\ 



(48) 
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The relation between a and E is exactly the same as in the Keplerian case, and a has a normal 
Keplerian meaning at a > 6. At r < b, however, the s-space a is very different; it can never be 
smaller than b and can be written as 

a = b^ > b. (49) 



Using dt = dr/v r , one finds 

At = 

2|£| V(s2-s)(s- Sl ) 



dt = J- - (s " 1)dg . (50) 



Let us define angle tp by 



S2 ~ si 

Then the time of motion from a given r to the pericenter r\ is 



sin tp = . (51) 



7T 2 

where T is the period of radial oscillations, 



T ( sin 2ib \ 
t P (r) = -U-e— ^ , (52) 



v^l^l 3 / 2 ' 

and e is the "s-space eccentricity" defined by 

_ fr(fi2 ~ si) 
6 ~ 2a 



T = 2t(r 1 ^r 2 ) = ^ n ,, (53) 



(54) 



The relation between E and T is the same as in the Keplerian case, and the formula for t p is 
formally the same (although the eccentricity e is defined in the s-space here). 

The eccentricity can also be expressed in terms of E and I, 

E\ 2 l 2 E 



^V-fo) + wei- (55) 

Linear orbits with I = have a maximum eccentricity 

e max (£) = l- J^. (56) 

In the limit r <^b (which corresponds to E = E — Eq <C \Eq\) the orbital motion occurs inside 
the homogeneous spherically-symmetric halo with density 
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